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Abstract 

Image restoration can be implemented efficiently by calculating the convolution of the digital 
image and a small kernel during image acquisition. Processing the image in the focal-plane in this 
way requires less computation than traditional Fourier-transform-based techniques such as the 
Wiener filter and constrained least-squares filter. In this paper, the values of the convolution 
kernel that yield the restoration with minimum expected mean-square error are determined using a 
frequency analysis of the end-to-end imaging system. This development accounts for constraints on 
the size and shape of the spatial kernel and all the components of the imaging system. Simulation 
results indicate the technique is effective and efficient. 


1 Introduction 

The Wiener filter is probably the best known and most widely used restoration tool. Given a few 
assumptions and some knowledge of the system, the Wiener filter minimizes the expected 
mean-square-error (MSE) of the restoration. While MSE is by no means a perfect yardstick for 
restoration quality, it is a. useful measure and leads to an optimal filter. In many applications, such as 
those requiring television-rate processing (30 images per second), the most serious drawback of the 
Wiener filter is its high computational cost. Small spatial kernels can be applied with much less 
computation. This paper describes the design of small restoration kernels that, within the spatial 
constraints, minimize restoration MSE. 


2 End-to-End Analysis and Wiener Restoration 

Traditionally, Wiener restoration has been based on a model of the imaging process with two 
components: the linear, shift- invariant point-spread function (PSF) of the image acquisition device 
and additive, signal-independent noise. This model ignores the significant impact of sampling and 
display reconstruction on image quality. A recent paper[l] presented a derivation of the Wiener filter 
that is based on a more accurate model of the end-to-end imaging process. This model is illustrated 
in Figure 1. 

The end-to-end process is described equivalently by equations in either the spatial domain or 
frequency domain. The displayed (or resulting) image r is 

1 i / 1 r oo \ 

r ( x ) = yj s(n- x')h(x')dx'N(n) + e[n]J f[n' - n]d(x - n') (1) 
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Figure 1: End-to-End Imaging and Spatial Restoration Model 


Assuming the scene s is periodic, the equivalent frequency domain expression for the spectrum of the 
result f is 


r v - 


E s[v']h[v']k[v - v'] + e[v\ f[v]d[v] 


( 2 ) 


where the notation f[u\ indicates the spatial frequency v/N, v cycles per N spatial units, of the 
Fourier transform of the image r. 

The Wiener filter minimizes the expected mean-square difference between the scene s and the 
resulting image r: 


S 2 = E ljjJ 0 K*) - r(z)| 2 cf 3 :j 

= e\ E MM-rMl 2 ) 


( 3 ) 


If the scene s and noise e are uncorrelated, stationary processes with power spectra 4> s and 4> e 
respectively, the expected mean-square restoration error can be rewritten in a form that is suitable 
for minimization: 

CO / CO 

s 2 = E ^ E - "'] 
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( 4 ) 


Minimizing the mean-square error with respect to the filter transfer function values f[v] leads to the 
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FIDELITY LOSS 


Reconstruction Blurs for DTF Factor = 0.3, 0.4, 0.5, 0.6 



Figure 12: Reconstruction Blur vs CTF FACTOR (Cat) 
(for DTF FACTOR = 0.3 to 0.6) 


Sampling Degradations for DTF FACTOR = 0.3, 0.4, 0.5, 0.6 



Figure 13: Sampling Blur vs CTF FACTOR (Cat) 
(for DTF FACTOR = 0.3 to 0.6) 


Image Blur and SR Degradation for DTF FACTOR = 0.3, 0.4. 0.5, 0.6 



Figure 14: Image and SR Blur vs CTF FACTOR (Cat) 
(for DTF FACTOR = 0.3 to 0.6) 


Fig 15 : Original Dollar Image with no image 
and SR degradations. 



Fig 16 : Reconstructed Scene with CTF FACTOR=0.3 
and DTF FACTOR=0.3 


Fig 17 : Reconstructed Scene with CTF FACTOR=0.3 
and DTF FACTOR=0.7 


22 


ORIGINAL PAGE IS 
OF POOR QUALITY 



definition of the optimal filter: 
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e $«iy] h[^'] m[i/ - 1/] + $ e M | e d[v'\ w[v - u'\ 


( 5 ) 


il/ — — OO 


This is the optimal digital filter given in Equation 26 of [l]. 

The mathematics of the following section is simplified by rewriting the expression for mean-square 
error in Equation 4 as 


S 2 = £ N [c[u] - b[v]f m [v] - bT[u]f[v] + a[v] \f[i 


where 


a[v] 

b[u] 
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Then, the optimal filter transfer function / given by Equation 5 is written: 

h[u\ 


f M = 


a[p] 


( 6 ) 


( 7 ) 

( 8 ) 
( 9 ) 


(10) 


3 Imposing Spatial Constraints 

In the derivation of the previous section, the Wiener filter is determined by an equation in the 
frequency domain: 

= b W\ (ii) 

The spatial equivalent of this frequency domain product is the spatial convolution: 

Lj2a[n-n']f[n'} = b[n] (12) 

n' 

where 

a[n] 
b[n\ 


1/ 

(13) 

E*MWT 

(14) 
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This convolution is equivalently expressed as a linear system of A equations in A variables (the 
spatial filter values). The system of equations can be expressed in matrix form: 

Af = b (15) 

where the N x N coefficient matrix A is 

A [n',n"} = Jj a [n' - n"] (16) 

the A X 1 result matrix b is the array b defined in Equation 14, and f is the A X 1 matrix of digital 
restoration PSF values to be determined. 

In the system of equations for the Wiener filter, there are as many equations as pixels in the image. 
However, if the size of the spatial restoration kernel is constrained, the system of independent 
equations can only be as large as the number of nonzero elements in the spatial kernel. The 
spatially constrained kernel is designed by specifying the system of linear equations whose solution 
will minimize mean-square-error within the constraints. 

The spatial constraint is expressed as a nonempty set of spatial locations, C, for which the 
restoration kernel can be nonzero. The elements that are not in the constraint set must be zero: 


f[n] = 0 if n <?C C {0... A - 1} 


(17) 


If all of the points in the restoration kernel are allowed to be nonzero (i.e., C = {0 ... A — 1}), then 
the optimal spatial kernel is the inverse transform of the Wiener filter (i.e., the solution of 
Equation 12 or 15). 

The expression for MSE is defined in Equation 6 in terms of the transfer function of the optimal 
filter. Before this expression can be minimized with respect to the restoration kernel values, it must 
be expressed in terms of those elements. The filter transfer function expressed in terms of the 
spatially constrained kernel values is 

/M = W/[«]»'r (is) 

nec 


Substituting this expression into Equation 6, yields the MSE in terms of the constrained kernel values: 
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(19) 
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( 20 ) 


Minimizing with respect to the restoration kernel elements yields 

TF J2 f[n']a[n - n'} = 6[».] n G C 

* n'eC 

This is an equation with a number of unknowns equal to the number of non-zero kernel values — \C\. 
There are \C\ equations (differentiating with respect to each of the constrained kernel elements) in \C\ 
unknowns (the |C| kernel values). This system of equations can be written as the matrix equation: 

Acfc = b c (21) 

where Ac is the \C\ X \C\ coefficient matrix, ic is the |C| X 1 matrix of kernel values, and be is the 
\C\ X 1 result matrix. 

The output matrix be of Equation 21 for the constrained filter is a submatrix of the corresponding 
matrix b of Equation 15 for the Wiener (unconstrained) fdter. The elements of the matrix be are the 
elements of b that are in the constraint set C . Similarly, Ac is a principal submatrix [2] of the 
coefficient matrix A consisting only of the rows and columns of A named in the constraint set C . 


4 Simulation Results 


This section presents restoration results for artificial scenes degraded by simulated imaging devices 
(as described in [3]). The problem design included two variables: the width of the acquisition transfer 
function and the noise level. Three cases for each variable were considered, producing a total of nine 
experimental restoration problems. Each of the nine problems was restored with kernels constrained 
to a number of sizes. Then, the accuracy of the constrained restorations was compared to the 
accuracy of the unrestored display and Wiener restoration. 

One-dimensional Fourier scenes were generated by specifying the spectral magnitude of a finite 
Fourier series and randomizing phase. The scene spectral magnitude s p was set to 




A' exp (|i/| /a s f aS ) if 0 < |;/| < 2 N 
0 otherwise 


( 22 ) 


with a s = N / 16 and (3 S = 0.75. Because the spectral magnitude is zero at the origin (s p [0] = 0), the 
resulting ensemble of scenes is zero-mean. The constant A' was defined as 0.0704946 so that the 
scenes had unit root-mean-square RMS energy. 

The model of the acquisition device transfer functions was suggested by Johnson[4]: 

h p [u] = exp (— (M /a h ) 0h ) (23) 

All three of the transfer functions in this section are bell curves (/3h = 2). 


• With a >h = 0.75, the transfer function roll off is mostly above the Nyquist limit. This function 
attenuates frequency components within the Nyquist limit only slightly and will therefore cause 
little blurring. However, the transfer function significantly passes components above the 
Nyquist limit and is therefore vulnerable to aliasing. 

• With Q/j = 0.50, the transfer function rolls off at a lower frequency and therefore causes 
somewhat, more blurring, but is less vulnerable to aliasing. 
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• With ah = 0.25, the transfer function is nearly zero beyond the Nyquist limit. This function 
virtually eliminates aliasing, but the resulting images may be blurred substantially. 

Three levels of zero-mean white noise were considered. Signal-to-noise ratio (SNR) is the ratio of 
RMS energy of the scene to RMS energy of the noise: 



For the low-noise images (high SNR), SNR— 100. For the moderate-noise images, SNR=25. For the 
high-noise images (low SNR), SNR=5. 

Real display devices are a significant component of the end-to-end imaging process but are not 
usually a source of much variability. Therefore, the simulated display function was not varied in these 
experiments — a single display model was used for all of the simulations. Schade[5] suggested a display 
model consisting of the sum of two Gaussian spots — the nucleus , a strongly-peaked central spot that 
contains most of the energy, and a broad flare spot around the nucleus. The composite display 
transfer function is 

d[y] = D\ exp (- (\u\ /q'i) 2 ) + £>2 exp (- (\u\ /a 2 ) 2 ) (25) 

The parameters for the functions are taken from Schade’s results: for the nucleus, D\ = 0.76 and 
Q'i = 0.4301484; for the flare, D 2 = 0.24 and a 2 = 0.0323814. For practical reasons, the display 
transfer function is cut off at twice the sampling rate ±2iV (the same length as the Fourier series used 
to generate the artificial scenes). The effect of the truncation is insignificant. 

Figure 2 illustrates the end-to-end imaging simulation for a representative scene. The top graph is 
the scene. Directly below it is the image created by applying the acquisition function with medium 
blur (ah = 0.50) to the scene. The third graph is the sampled image. Next is the sampled scene plus 
moderate noise (SNR = 25). The bottom graph of Figure 2 shows the unrestored display. Acquisition 
blurring, aliasing due to sampling, additive sensor noise, and display degradation are all present in the 
output of the system. The goal of restoration is to process the noisy digital image shown in the fourth 
graph so that when it is displayed, the output (the bottom line) is more like the input (the top line). 

The spatial kernels were constrained to have zero value at all but an odd number of locations 
centered at the origin — the smallest kernel, with three elements, was allowed non-zero values only 
where |?t| < 1; the next smallest, with five elements, was allowed non-zero values only where |n| < 2; 
and so on. The largest constrained kernel has (N - 1) elements; only the element at n = N/2 was 
constrained to 0. The next-largest optimal kernel (no elements constrained to 0) is the spatial kernel 
of the Wiener fdter. 

The optimal three-point and five-point, kernels for the example of Figure 2 and the corresponding 
transfer functions are shown in Figure 3. The Wiener filter transfer function and part of the 
corresponding spatial kernel are also illustrated. Only the first few elements of the Wiener kernel are 
shown; the magnitude of the Wiener kernel elements beyond 6 pixels from the origin is less than 
O.OliV. Clearly, the optimal small kernels are quite different than the kernels produced by a 
truncating the Wiener PSF. As can be seen by comparing the transfer functions, the optimal 
three-point kernel does a. fair job of approximating the Wiener filter at low frequencies but amplifies 
high-frequency components where SNR is lower much more than does the Wiener filter. The transfer 
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function of the optimal five-point kernel more closely approximates the Wiener filter, but is still quite 
different. 

Figure 4 shows the original scene, the unrestored output, the output with three-point restoration, the 
output with five-point restoration, and the output with Wiener restoration. Visual comparison is a 
subjective process, but it is clear that all of the restorations are more like the original scene than the 
unrestored output. It is more difficult to conclude from visual inspection which of the restorations is 
the best. Some of the features seem to be restored best by the three-point kernel; other features are 
best restored by the Wiener restoration. 

Figure 5 presents numeric measures of restoration accuracy as a function of kernel size. Restoration 
accuracy is described by the RMS difference between the displayed image and the scene, relative to 
the RMS energy of the scene: 

E S - '’HI 2 

Relative RMS Error = P — (26) 

^ 2_. I S MI 

Each of the nine restoration problems was performed 32 times — that is, each execution used a 
different scene from the ensemble and different random noise. The plots show the relative RMS error 
averaged over all 32 executions. The standard deviations of the relative RMS error were so small that 
plotting them on these graphs proved impractical. The plots are shown only for kernels with 65 
elements or fewer (radius 32). In all cases, only negligible improvement occurred beyond 19 elements 
(radius 9). (The kernel of the Wiener filter lias 255 elements, a radius of 128.) 

In many cases, the three-point and five-point kernels yielded results that are nearly as accurate as the 
Wiener filter. This is particularly true when there is little noise (e.g., SNR=100 — the leftmost 
column). Small kernels are relatively less successful in low SNR situations (e.g., SNR=5 — the 
rightmost column). In low SNR problems, the restoration kernel should suppress noise by local 
averaging, but small kernels are restricted in doing so by their size. 

In the image with medium blur and medium noise (ah = 0.50 and SNR=25) — the middle graph of 
Figure 5 — the average unrestored relative RMS error was 0.204613. The Wiener filter reduced this to 
0.051149, a decrease of 0.1536464 or 75%. The three-point kernel resulted in an error of 0.091685, a 
decrease of 0.112928 or 55%. The three-point kernel (radius 1) achieved 73% of the improvement of 
the Wiener filter. The five-point kernel (radius 2) reduced the relative RMS error to 0.083614, a 
decrease of 0.120999 or 59%. This is 79% of the improvement of the Wiener filter. These small 
kernels achieve a large portion of the improvement of the Wiener filter with less computation. 


5 Conclusion 

Restoration implemented by convolution with a small kernel requires less processing than traditional 
Fouricr-transform-based techniques such as the Wiener filter. Because convolution with a small kernel 
is a local operation, it is easily applied in parallel on all pixels in the focal-plane during image 
acquisition. The simulation results indicate that the optimal constrained restoration kernel effectively 
restores continuous, one-dimensional functions degraded by blurring, sampling, noise, and 
reconstruction — the types of degradations found in real imaging systems. Similar results were 
observed in simulations not presented in this paper using other one-dimensional scenes with different 
statistics. Two-dimensional simulations and actual restorations are in preparation. 
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